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Phase diagrams in the three-flavor Nambu-Jona-Lasinio model with the Polyakov loop 
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We present extensive studies on hot and dense quark matter with two light and one heavy flavors 
in the Nambu-Jona-Lasinio model with the Polyakov loop (so-called PN JL model) . First we discuss 
prescription dependence in choosing the Polyakov loop effective potential and propose a simple and 
rather sensible ansatz. We look over quantitative comparison to the lattice measurement to confirm 
that the model captures thermodynamic properties correctly. We then analyze the phase structure 
by changing the temperature, quark chemical potential, quark masses, and coupling constants. 
We particularly investigate how the effective Ua(1) restoration and the induced vector-channel 
interaction at finite density would affect the QCD critical point. 

PACS numbers: 12.38.Aw, 11.10. Wx, 11.30.Rd, 12.38.Gc 



INTRODUCTION 



The phase diagram of hot and dense matter out of 
quarks and gluons described by Quantum Chromody- 
namics (QCD) has attracted theoretical and experimen- 
tal interest for decades P, H| . We can define one phase 
transition associated with chiral symmetry restoration 
in the vanishing quark mass limit (i.e. m q — > 0) which 
is commonly referred to as the chiral phase transition. 
In the quenched limit with infinitely heavy quark mass 
(i.e. m q — > oo), on the other hand, we can define an- 
other phase transition from the hadron (glueball) phase 
to the color deconfinement phase. The question is then 
the nature of these phase transitions with intermediate 
quark masses of two light (up and down) and one heavy 
(strange) flavors. We stress that the chiral and decon- 
finement phase transitions are conceptually distinct phe- 
nomena and, theoretically speaking, they reside in the 
opposite limits with respect to the quark mass. Never- 
theless, the standard QCD phase diagram on the plane 
of the temperature T and the chemical potential /i has 
only a single transition or crossover boundary. Whether 
this is really the case is not trivial a priori and not quite 
settled yet. 

It is the result from the Monte-Carlo integration of the 
(quenched) QCD partition function on the lattice that 
had led us to this phase diagram with a single phase 
boundary Q. (See Ref. [H and references therein for 
historical background.) Later on, the lattice QCD sim- 
ulation with dynamical quarks [5[ have confirmed that 
the chiral and deconfinement phase transitions occur at 
the same temperature (or at different but close tempera- 
tures @). This observation suggests that two phenomena 
of chiral restoration and color deconfinement should be 
locked by some dynamical mechanism so that they should 
take place (nearly) at once. 

We can find the first successful study based on dynam- 
ical model to give an account for this locking mechanism 
in the work by Gocksch and Ogilvie [7j. They have con- 
structed the effective action of QCD by means of the 
strong coupling and large dimensional expansions. The 
same action has been discussed at finite T and \x also 



by Ilgenfritz and Kripfgantz [8]. There were proposed 
some generic mixing arguments which aim to explain the 
coincidence of critical temperatures in a model indepen- 
dent way @. The mixing effect, however, does not suf- 
fice to force two crossovers to be a single one in view 
of the associated peaks in the susceptibility; two sep- 
arate crossovers (susceptibility peaks) with mixing for 
each cannot be ruled out. That means, the mixing ef- 
fect is a necessary but not sufficient condition in order to 
realize the coincidence in a way seen on the lattice [Icj . 
Thus, the locking between chiral restoration and decon- 
finement should need more tangled dynamical properties 
of two phenomena. 

To reveal the relevant dynamics, the present author 
proposed a useful model II 311 based on the Nambu-Jona- 
Lasinio (NJL) model [TTI. Il2l| with the Polyakov loop de- 
grees of freedom augmented, which was inspired by the 
strong coupling analyses @, 0, [H| ■ The peculiar feature 
in this model is that we can uniquely determine the cou- 
pling between the chiral condensate, which is an order 
parameter for the chiral phase transition in m q — > 0, and 
the Polyakov loop, which is an order parameter for the 
deconfinement phase transition in m q — > oo. The model 
inputs and outputs have been carefully compared to the 
lattice QCD data by Ratti, Thaler, and Weise [l5[ and 
they named this hybrid description as the PNJL model. 

The PNJL model has been generalized to the three- 
flavor case recently [H, Qj} • In the present paper we shall 
extensively explore the phase diagrams in the three-flavor 
PNJL model by changing four physical variables, namely, 
T, \i, the light quark mass m u( i, and the heavy quark 
mass m 8 . First we revisit the choice of the Polyakov 
loop effective potential that cannot avoid ambiguity in 
the PNJL model approach. We claim that a careful con- 
sideration is necessary for the effective potential form. 
Once we fix the pure gluonic sector by specifying the 
potential, we can calculate the mean-fields of the chiral 
condensate and the Polyakov loop to draw the phase di- 
agrams. 

Although it is usually assumed implicitly, we have no 
first-principle insight into the locking of two crossovers in 
the finite-density region. The lattice Monte-Carlo simu- 
lation is of no practical use except when fi is much smaller 
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than T. So far it seems that almost nothing but the PNJL 
model can access both transitions at any density. Strictly 
speaking, in fact, the mean-field treatment of the PNJL 
model is not totally free from the sign problem. Detailed 
analyses in Ref. [l8| support that the saddle-point of the 
mean-field energy leads to an appropriate estimate for 
the mean-fields, however. Hence, we will not argue the 
sign problem any more in the present work. 

In this paper, after we check that our results from the 
three-flavor PNJL model are consistent with the state-of- 
art lattice simulation at zero chemical potential [l9l |, we 
will shift our emphasis toward the QCD (chiral) critical 
end-point in the last half. It should be noted that the 
terminal point of the first-order phase boundary has a 
second-order phase transition characterized by the uni- 
versality class of the Z(2) Ising model and this special 
point is often called the QCD critical point. The search 
for the critical point is one of the most interesting prob- 
lems in finite density QCD because it provides us with 
a firm milestone for our quest for the QCD phase dia- 
gram. If we are lucky enough to find out the critical 
point as predicted in theory, we can get confident about 
our theory reliability. This is, so to speak, a mutual cor- 
respondence between theory and experiment, which is an 
ideal situation for sound scientific developments. 

The existence of the QCD critical point has not been 
established yet. We cannot exclude a possibility that the 
QCD phase transition is smooth everywhere in the fi-T 
plane, while there are a pile of indirect evidences for its 
existence [H [H [H [H IM [H, [1| . In model studies, in 
fact, a minor modification in the treatment could easily 
smear a first-order phase transition out into a crossover, 
as demonstrated later. In particular we shall pay at- 
tention to two obscure factors which may significantly 
affect where the critical point is and even whether it ex- 
ists. Those two factors are the magnitude of the Ua(1)- 
breaking anomaly interaction and the vector-channel in- 
teraction. The former, the UA(l)-brcaking term, in- 
duces a six-quark vertex called the 't Hooft term which 
mixes three different flavors up and is responsible for the 
first-order phase transition in the chiral limit [271 ]. It 
could be possible at finite temperature and density that 
the 't Hooft interaction is reduced by instanton suppres- 
sion [23, [23, Ho] . The latter effect, i.e. the vector-channel 
interaction term, does not break chiral symmetry and the 
zcroth component directly couples the quark density. It is 
thus likely that the finite-density environment enhances 
or induces interactions in the vector channel which could 
weaken the first-order phase transition [2(| [3l|, [H, [33| ■ In 
this paper we shall quantify these effects on the location 
of the QCD critical point using the three-flavor PNJL 
model. 



II. MODEL SETUP 

The present author proposed the PNJL model action 
in Ref. [l3| inspired by the effective action in strong- 



coupling QCD with dynamical quarks 0, H, It is 

possible to some extent to elaborate a field-theoretical 
setup for the PNJL model starting with the Lagrangian 
density [TH ]. For this purpose it is required to assume 
a homogeneous mean-field distribution of the Polyakov 
loop. In other words, the temporal component of the 
gauge field, A4, in Euclidean space-time must be approx- 
imated by a spatially constant mean-field, so that one can 
perform the one-loop integration with respect to ther- 
mal quarks in a closed form. This thermal integration 
leads to the unique coupling between the chiral conden- 
sate and the Polyakov loop. Spatial uniformity is in fact 
a mean-field ansatz, however, and it makes a contrast to 
the strong-coupling framework as we shall discuss 
shortly. 

In the PNJL model the Polyakov loop is therefore put 
in as a global mean-field rather than a local dynamical 
variable, which is analogous to the treatment of the chiral 
condensate in the ordinary NJL model; the Lagrangian 
density with a shift by the mean-field is sometimes re- 
ferred to as the mean-field Lagrangian that contains no 
kinetic term for the mean-field. Such an approximation 
should work to investigate the bulk property of the ther- 
modynamic system, while we have to be aware that the 
mean-field model cannot properly deal with the spatial 
structure of confined objects. It is beyond the scope of 
the simple PNJL model framework, for instance, to ex- 
tract the heavy-quark potential. 

All the model ingredients are thus given as mean-field 
variables. Here we would prefer to start with the mean- 
field free-energy after one-loop integration for the model 
setup. Let us decompose the free-energy below into four 
pieces and discuss them in order. That is, the total free- 
energy (or the grand potential) is a sum of four contri- 
butions; 

^PNJL = ^cond + ^ zero ^quark ~l~ ^Polyakov , (1) 
NJL part pure gluonic part 

where ri C ond represents the condensation energy in the 
chiral sector, Q zcro the zero-point energy which is impor- 
tant in the NJL model formulation, f2 qua rk the thermal 
quark contribution with the Polyakov loop coupling com- 
ing from the Dirac determinant, and Opoiyakov gives the 
effective potential in terms of the Polyakov loop variable. 
As indicated in Eq. (TTJ , we can deduce the first three from 
the standard NJL model and the last one from the pure 
gluonic theory. 

A. Condensation Energy 

We can read the condensation energy from the stan- 
dard NJL model Lagrangian. Using the notation by Hat- 
suda and Kunihiro fl2| , we write the four-quark inter- 
action in the scalar channel and the six-quark 't Hooft 
interaction as 

£s = f [(^) 2 + feAaV) 2 ] , (2) 
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and 



C. Thermal Quark Energy 



Ca = 3D [det lp(l - 75)^ + h.c] , 



(3) 



respectively. For later convenience we also give an ex- 
pression for the vector-channel interaction; 



(4) 



For the moment we will work only in the gy = case. 
Here, A a 's ar e th e Gell-Mann matrices in flavor space 
(with Xq = y/2/3) and the matrix determinant is taken 
also in flavor space. In the mean-field approximation 
with three condensates, (in), (dd), and (ss), the scalar 
four-quark interaction is rewritten as 

gs(uu) 2 — > gs(uu — (uu) + (uu)) (uu — (uu) + (uu)) 
~ gs(uu) 2 + 2gs(uu)(uu — (uu)) 
= -g s (uu) 2 + 2g s (uu)uu, (5) 

in the w-quark sector and likewise for d-quarks and s- 
quarks. In this way we can readily reach the following 
expression for the condensation energy; 

ftcond = gs((uu) 2 +{dd) 2 +(ss) 2 )+Ag D {uu}{dd){ss) . (6) 

We see that the six-quark interaction induces the flavor- 
mixing interaction indeed which makes the phase transi- 
tion of first-order in the presence of massless three flavors. 



B. Zero-Point Energy 

The zero-point energy diverges and requires the ultra- 
violet cutoff A to regularize the three-momentum inte- 
gration. Since the NJL model is a non-renormalizable 
cutoff theory depending on the choice of A, the zero- 
point energy contribution largely affects the model out- 
put. With the quasi-quark energy dispersion relation, 
£i(p) = \/v 2 + Mf , the zero-point energy can be ex- 
pressed simply as a summation of all £j(p)/2, that is, 



d 3 p 
(2tt) 3 



£i(p) , 



(7) 



where 2 is the spin degrees of freedom, N c = 3 is the 
number of colors, and the particle and anti-particle con- 
tributions cancel 2 in the denominator of £j(p)/2. The 
constituent quark mass is defined as a sum of the current 
quark mass and the mean-field as 

M u = m u - 2g s (uu) - 2g- D (dd)(ss) , 

M d = m d - 2g s (dd) - 2g B {ss) (uu) , (8) 

M a = m s - 2g s (ss) - 2gu{uu)(dd) , 

which is understood from the second term in Eq. l[5]). 



The thermal quark energy is where we can uniquely in- 
troduce coupling between the chiral condensate and the 
Polyakov loop. In the PNJL model, under the assump- 
tion of the presence of the spatially uniform Polyakov 
loop background, the one-loop free-energy is modified as 



n, 



quark 



2T£ 



d 3 P 



{ 



In det 



(2tt) 3 

lndet[l + it e -fe(rt+M)/T]| 



(9) 



Let us comment on preceding works [3J, l35l | in which a 
similar coupling form is addressed. 

We note that the above expression is identical with 
that in the strong coupling expansion but the physics 
content is slightly different. The Polyakov loop L is a 
mean-field from the beginning here, whereas the strong 
coupling calculation at finite temperature decouples the 
temporal hopping from spatial link variables [141] . As 
a result, the quark excitation is static in the strong- 
coupling leading order, and the above expression results 
at each lattice site in this way, that is, L could be a local 
variable in the strong coupling expansion. 

It is noteworthy that the three-momentum integration 
above is finite and has no need for the ultraviolet cutoff. 
We can thus relax the cutoff in the thermal quark energy, 
though we found that the s-quark sector behaves unnat- 
urally at extremely high temperature without the cutoff, 
which is of no importance practically. In this work we 
will not impose the momentum cutoff onto the thermal 
quark energy in order to let the thermodynamic quanti- 
ties free from cutoff artifact. 

The Polyakov loop L is an N c x N c matrix in color space 
and is defined originally in terms of A4. The explicit form 
of the Polyakov loop is irrelevant in our study because 
we treat it as a model variable and will not return to the 
original definition of the Polyakov loop in terms of the 
gauge field. 

In the simplest mean-field approximation one can ex- 
press the free-energy as a function of the traced Polyakov 
loop expectation value defined by 

£=—(trL), i=—(txL^). (10) 

It should be mentioned that we must distinguish I and 
I at finite density [H, [H, |37|; both I and i are real, 
and nevertheless, I > I whenever \i > 0. This is 
because a finite chemical potential gives rise to a C- 
odd term like u ImftrL] in the average weight leading 
to I - I - ^((Im[trL]) 2 ) > for_ small u. We can 
also give an intuitive explanation; I represents the ex- 
ponential of the free-energy gain, fi, by the presence of 
an anti-quark. The test charge brought in by an anti- 
quark can be more easily screened in a medium with more 
quarks than anti-quarks. Therefore, fj < fi, that means, 
I = b~M t > I = e~^/ T for a positive /x. 
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It is straightforward to take an average of the 3x3 
determinant to reach 



det [l + L c-( £ -" )/ T ] )=l + e-^-^' T 

+ 3£e -(e-^/T + 3le -2(e- tl )/T^ (n) 

det [1 + L* e -(^)/T] } = 1 + c -3( £ +m)/t 

+ 3^ e -(e+p)/T + 3£e -2( £+A1 )/T_ (12) 



In this work we use the logarithm of the above expressions 
as the mean-field free-energy and will not perform the 
group integration over L. Roughly speaking, the approx- 
imation involving the group integration [7| corresponds 
to what is called the Weiss mean-field approximation in 
the spin system. The integration has an effect on the 
quantitative results [38| but a simple mean-field treat- 
ment suffices for our present purpose. We also remark 
that the action is invariant under simultaneous replace- 
ment I <-> I and — fi <-» 



D. Polyakov Loop Energy 

In the definition of the PNJL model the choice of the 
Polyakov loop potential has subtlety because the effective 
potential has not been known directly from the lattice 
QCD simulation. In the present study we will assume 
the strong-coupling inspired form of 



^Polyakov = -6-r{54e" a / T II 

+ ln[l-6«-3(«) 2 + 4(t 3 + £ 3 )]} 



(13) 



The logarithmic term appears from the Haar measure of 
the group integration with respect to the SU(3) Polyakov 
loop matrix. The first term is reminiscent of the near- 
est neighbor interaction in the effective action at strong 
coupling. The temperature-dependent coefficient of this 
11 term controls the deconfinement phase transition tem- 
perature. It should be, however, noted that the model 
parameters are assumed to be temperature- independent. 
(See Ref. [3^| for the running coupling effects including 
renormalization.) 

In this simple ansatz for the Polyakov loop potential, 
we have two parameters; a and b. The deconfinement 
phase transition is determined solely by the choice of a, 
while b parametrizes the relative strength of mixing be- 
tween the chiral and deconfinement phase transitions. If 
b is small, chiral restoration dominates the phase transi- 
tion, and if b is large, deconfinement is more governing. 

We will numerically make a comparison between the 
above-proposed ansatz and others in the next section. 



III. NUMERICAL PROCEDURES 

Now that we have specified all the constituents in the 
model action, we get ready to proceed to the numerical 



analyses. We will solve the following four equations in a 
self-consistent way, 



dfl 



PNJL 



PNJL 



on 



PNJL 



an 



PNJL 



d{uu) d{ss) 



Of 



Of 



(14) 



to acquire (uu) = (dd), (ss), t, and I as functions of 
the model input. For this purpose we have to fix all the 
model parameters, A, gs, 3d in the NJL potential, and a 
and b in the Polyakov loop potential. 



A. Parameter Choice 

The Polyakov loop coupling appears only in the ther- 
mal part, that means that the NJL model parameters 
fixed at T = ji = are not affected by introduction of 
the Polyakov loop coupling. In this work we will employ 
the widely accepted parameter set according to Hatsuda 
and Kunihiro [12j |: 



A = 631.4 MeV, 
mud = 5.5 MeV, 
5S -A 2 = 3.67, 



m s = 135.7 MeV, 
g D ■ A 5 = -9.29 , 



(15) 



which nicely reproduces the tt mass, the K mass, the rj 
mass, and the ir decay constant f v . Here m u d stands rep- 
resentatively for the light quark mass, i.e. m u d = rn u — 
m d . 

Regarding the Polyakov loop potential, we can fix the 
parameter a by the condition that the first-order phase 
transition in the pure gluodynamics takes place at T = 
270 MeV, which yields 



a = 664 MeV , 



(16) 



and then the remaining variable is b only. Actually, the 
determination of b suffers uncertainty and there is no 
established prescription. In this study we shall take a 
value of b that leads to simultaneous crossovers of chiral 
restoration and deconfinement around T ~ 200 MeV. As 
a result, we set 



b ■ A" 



0.03. 



(17) 



B. Other Polyakov Loop Potentials 

The choice of the Polyakov loop potential has some 
variations, as we have mentioned before. Our choice of 
Eq. (fT3"|) is much simpler than the widely accepted forms 
by Ratti, Thaler, and Weise [l5| and by Rofiner, Ratti, 
and Weise [40l |. It would be instructive to scrutinize re- 
spective forms and quantify the difference numerically. 
Let us call the "RTW05 potential" to indicate 



RTW05 



(18) 
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with b 2 (T) =a + ai (To/TV + a 2 (T /T) 2 + a 3 (T /T) 3 , 
which is proposed in Ref. [15| . There are seven param- 
eters, a Q — 6.76, a\ = —1.95, a 2 = 2.625, a 3 = —7.44, 
6 3 = 0.75, 64 = 7.5, and T = 270 MeV such that the 
potential (| 18[) reproduces the pressure, energy density, 
and entropy density in the pure gluonic sector measured 
on the lattice. A slightly different choice is suggested in 
Ref. [H which we shall call the "RRW06 potential" ; 

^RRW06 = T A { 

1 , (19) 

+ b(T) ln[l - 6U - 3{££) 2 + 4(£ 3 + P)] I 

with a(T) = a + ffli(T /T) + a 2 (T Q /T) 2 and b(T) = 
&3(To/T) . Five parameters are fixed as ag = 3.51, at = 
-2.47, a 2 = 15.2, b 3 = -1.75, and T = 270 MeV. We 
note that 63 plays the same role as b in our ansatz (fl~3|) . 
Actually, if we substitute T = 190 MeV to lower the 
crossover temperature as argued in Ref. , b^T§ ■ A 3 ~ 
0.044 (where A is not our value but 650 MeV used 
in Ref. [l5j]) which turns out to be comparable to our 
choice (fTTj) . 

Under the assumption that fircrwos an d ^RRW06 cor- 
respond to the total negative pressure in the pure gluonic 
theory, they approach the Stefan-Boltzmann limit at high 
temperature, that is, p = (2 -8 • n 2 /90)T 4 = 1.75T 4 . One 
can easily make this sure from — ao/2 — 63/3 + 64/4 = 
-1.75 in S! RT wo5 and -ao/2 = -1.75 in S! RRW06 . 

We would claim, however, that f2RTW05 an d ^rrwo6 
might overcount the relevant degrees of freedom in the 
system. In the high temperature limit not only the 
Polyakov loop but also the deconfined transverse glu- 
ons contribute to the pressure. Since the Polyakov loop 
corresponds to the longitudinal gauge field, the Stefan- 
Boltzmann limit should be saturated by the transverse 
gluons but not the Polyakov loop. It is thus a subtle as- 
sumption that the effective potential with respect to the 
order parameter field can reproduce the total pressure, 
energy density, and entropy density for all temperatures. 

One can understand this from a more familiar exam- 
ple. Let us consider the mean-field effective potential 
in the 0(4) linear sigma model. The effective poten- 
tial with respect to the a condensate describes the chiral 
phase transition. The total pressure should contain con- 
tributions from the ir excitations too which are not fully 
included in the effective potential in terms of (a). 

It is not our point to insist that ^rtwos and ^rrwo6 
are doubtful. Our main point lies in the other way around 
in fact. We presume that their parametrization works 
in effect for the following reason; the pressure contribu- 
tion from transverse gluons is a function of T, and the 
Polyakov loop is also a function of T, and so the former 
can be implicitly parametrized by the latter. Then, it 
is possible to express the total pressure in the form of 
Eq. (fl8|) or (fl~9|) . One has to keep in mind, however, 
that the total pressure in this interpretation would make 
sense provided that the Polyakov loop is already solved as 
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FIG. 1: Comparison of the Polyakov loop pressure excess as 
a function of the temperature. The vertical axis signifies the 
effective degrees of freedom. RTW05 and RRW06 represent 
Ortwos and fiRRW06 with £ and £ given as a solution of the 
full gap equation with two quark flavors. For comparison we 
used the same NJL model parameters, A, gs, and m u = mj 
as in Ref. [l5| and quenched the s-quark sector to draw the 
solid curve which represents our Opoiyakov 



a function of T. Therefore, one should solve Eq. (|14[) first 
and then one can fit the total pressure using Eq. (fl8|) or 
(fig)) with solved £(T) and £(T) substituted. One should 
not use the total pressure itself to optimize the varia- 
tional parameters £ and £. This may explain why the 
critical temperature determined with Eq. (fl8|) or (fl"9]) 
put into the gap equations becomes relatively higher. 
The Polyakov loop effective potential which overcounts 
the gluonic degrees of freedom would drag the crossover 
point closer to the pure gluonic transition temperature 
T = 270 MeV. 

We emphasize that our simple choice of the Polyakov 
loop potential is physically natural and, interestingly, 
it makes only little difference from the numerical re- 
sults based on the RTW05 or RRW06 potential. This 
sounds very good, for our new potential ansatz does not 
ruin the nice agreement to the lattice data addressed in 
Refs. HE S3, In Fig. [0 we plot the Polyakov loop 
pressure difference from the zero temperature value us- 
ing the mean-fields obtained from the full gap equations 
with two flavors. We could, of course, show the genuine 
total pressure with the Polyakov loop and quark contri- 
butions both. We subtracted the quark contribution in 
Fig. [1] because the quark contribution makes the com- 
parison blurred; the quasi-quark pressure is dominating 
up to near T c but it is not sensitive to the choice of the 
Polyakov loop with r2 qua rk untouched. The non-trivial 
part is the Polyakov loop contribution that we now focus 
on here. 

In the absence of interaction, the pressure is given by 
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the Stefan-Boltzmann law, 7r 2 T 4 /90, multiplied by the 
effective degrees of freedom which we denote by v. To 
see how v increases as T goes up, we normalize the pres- 
sure by the Stefan-Boltzmann unit; 7r 2 T 4 /90. Clearly 
both JIrtwos an d ^rrwo6 increase with increasing T and 
asymptotically approach the value of v — 2 (polarization) 
x 8 (color)=16. It is so by construction, as we explained. 
It is intriguing to note that our ansatz (|13p results in the 
solid curve in Fig. [1] which is close to the dashed and 
dot-dashed curves by ^rtwos an d ^rrwo6 as long as the 
temperature is below 300 MeV ~ 1.5T C . We do not have 
to care much about the discrepancy in the higher tem- 
perature region, in fact, because the validity region in the 
present study extends at best up to ~ 2T C above which 
transverse gluons should be dominant. Therefore, we can 
conclude that all these potential choices are consistent to 
each other within the validity range of the temperature. 
In our choice (fT?j) the effective degrees of freedom slowly 
decrease at higher temperature in the Stefan-Boltzmann 
unit. This is reasonable because the Polyakov loop must 
give way to transverse gluons. 

The nearly coincidence of three curves in the vicinity 
of T c in Fig. [1] delivers us an important message. The 
Polyakov loop takes on a major fraction of the system 
pressure up to the temperature around 1.5T C . We should 
recall that two parameters, a and b, in Eq. (|13p have been 
fixed not to reproduce the pressure but just to yield To = 
270 MeV in the pure gluonic sector and T c ~ 200 MeV 
with 2 + 1 flavors. 



IV. ZERO DENSITY RESULTS 

Here we show the model results at zero quark density 
with our choice of the model parameters. In our subse- 
quent discussions we will make clear the virtues of the 
PNJL model as well as some caveats. 



A. Order Parameters 

Because nothing breaks isospin symmetry in this work, 
we will show the numerical results only for the u-quark 
sector which is degenerate to the d-quark sector. 

First of all, we present Fig. [2] to confirm that simulta- 
neous crossovers of deconfinement and chiral restoration 
certainly realize in the PNJL model. The chiral con- 
densates are normalized by their zero-temperature value; 
(uu)o = (246 MeV) 3 and (ss} = (267 MeV) 3 for light 
and heavy quarks, respectively. 

The reason why we find the simultaneous crossovers 
around T c ~ 200 MeV (the temperature derivative gives 
T c = 204.8 MeV) is that we have chosen the value of b 
as Eq. (fTT|) to adjust the crossover temperature by hand. 
Thus, we note that the crossover temperature is not the 
model output but the input. Nevertheless we would com- 
ment on a non-trivial feature inherent in the model dy- 
namics; the chiral phase transition can never occur until 




T [MeV] 

FIG. 2: Order parameters at zero density as a function of 
the temperature. The chiral condensates are normalized by 
(twi)o = (246 MeV) 3 and (ss} = (267 MeV) 3 . The solid 
and dotted curves represent the it-quark and s-quark chiral 
condensates, respectively, and the dashed curve represents the 
traced Polyakov loop I. 

the Polyakov loop grows up [HI]. It is also interesting 
to look at the behavior of the s-quark chiral condensate 
depicted by the dotted curve. The results for (ss) are the 
output rather than the input unlike (uu). If we define a 
crossover temperature for the s-quark sector, it should 
be higher than the simultaneous crossovers due to the 
explicit breaking of chiral symmetry. 

B. Effective Confinement 

Let us elucidate how the effective confinement is pos- 
sible in the model description. The underlying idea in 
the PNJL model is that the group integration (average) 
with respect to the Polyakov loop acts as a projection 
onto the center symmetric state (or the canonical ensem- 
ble [42| with zero Z3 charge) if there is no Polyakov l oop 
mean-field. (See Eqs. ([HI) and (O and also Refs. [HHl 
for details.) We solve the four coupled equations (fT4"|) at 
T ^ and [i — 0, and plot the quark pressure differ- 
ence from the zero temperature value in Fig. [3] using the 
obtained mean-fields. 

In the limit of massless two and three flavors we can 
count the fermionic degrees of freedom as v = (7/8) ■ 
3 • 2 • 4 = 21 and v = (7/8) • 3 ■ 3 ■ 4 = 31.5, re- 
spectively. Because the system of our interest is quark 
matter with two light and one heavy flavors, v should 
take a certain value between 21 and 31.5 at temperature 
above T c where chiral symmetry is restored. This ex- 
pectation is manifest in view of Fig. [3] both in the NJL 
model and in the PNJL model. Here we have determined 
the pseudo-critical temperature by the location where the 
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FIG. 3: Effective degrees of freedom associated with 2 + 1 
flavor quarks with increasing temperature in the NJL model 
and in the PNJL model. The critical temperature is T c = 
1 71 .6 MeV in the NJL model case and T c = 204.8 MeV in the 
PNJL model case. 



temperature derivative, d(uu)/dT 1 is largest. It follows 
that T c = 171.6 MeV for the NJL model results and 
T c = 204.8 MeV for the PNJL model results (see also 

Fig. m. 

Even in the standard NJL model the effective degrees 
of freedom go down as the temperature becomes lower. 
This is because quark excitations are suppressed by the 
constituent quark mass in the low temperature side where 
chiral symmetry is spontaneously broken. In reality the 
system should be mainly composed of a gas of tt° and tt ± 
below T c but the it mass is ~ 135 MeV which is compa- 
rable to the critical temperature. It is thus expected that 
we can neglect the 7r loop corrections in the pressure in 
the first approximation. 

We can see from Fig. |3] that the NJL model contains 
too many (unphysical) quark excitations below T c . It 
should be mentioned that, strictly speaking, too many 
excitations cannot be concluded until this comparison 
and the observation that the PNJL model is consistent 
with the lattice results. These fictitious excitations di- 
minish only slowly. It is apparent that the Polyakov 
loop projection works efficiently in the PNJL model case. 
The effective degrees of freedom rapidly decrease near T c , 
that means that artificial quark excitations are removed 
by the Polyakov loop coupling. Therefore, we can an- 
ticipate that the PNJL model should be more capable to 
capture realistic thermodynamics than the standard NJL 
model especially at temperatures near T c . Also, because 
the Polyakov loop projection affects the quark sector, it 
is a natural expectation that the PNJL model would be 
a more suitable description than the NJL model in the 
finite density region where quarks exist abundantly. 

Finally we shall remark that the separation of the total 



FIG. 4: Chiral and Polyakov loop susceptibility at zero den- 
sity as a function of the temperature in the 2 + 1 flavor PNJL 
model. 

pressure into the Polyakov loop and the quark contribu- 
tions like in Figs. [1] and [3] does not make sense in the 
mean- field approximation. This is because each of (uu), 
(ss), £, and £ determined by the gap equations (|14[) in- 
volve entangled contributions and thus a clear separation 
is impossible in any way. 



C. Susceptibility 

In this subsection we clarify how we can evaluate the 
susceptibility in respective channels of our interest. A 
useful alternative is to be deduced from the temperature 
slope, i.e. —d(uu) /dT, di/dT, and so on. They behave as 
a function of T in a similar manner as the susceptibility 
but the temperature slope is not really informative more 
than the order parameter curves read from Fig. [5] 

In order to compute the susceptibility in the mean-field 
model we need some caution. Remembering that the 
logarithm of the partition function is — VQpnjl/T, we 
can give the definition of the dimcnsionless susceptibility 
of our interest as 



1 d 2 { 


—Qpnjl/T) 


(20) 


4T 


dm ld 


1 d 2 {~ 




(21) 


T 


dm 2 . 




drj dfj 


(22) 


1 d 2 {- 


^pn.tl/7 1 ) 


(23) 


T 





for the light-quark susceptibility, the heavy-quark sus- 
ceptibility, and the Polyakov loop susceptibility, respec- 
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tively. We also enumerate the quark number suscep- 
tibility that we will discuss later. Here we have in- 
serted the Polyakov loop source i] and fj in the poten- 
tial as S! Polyakov — * fipoiyakov - T(r]£ + fj£). It is cru- 
cial to notice that we have to take the derivative in a 
way that it hits the mean-fields also. That means that 
we should take d(uu)/dm u d, d 2 {uu) / dm^ d , d(ss)/dm u d, 
d 2 (ss)/dm 2 td , d£/dm u d, d 2 £/dm 2 d , etc into account to 
evaluate Eq. ((20)) for instance. Otherwise we would miss 
the loop effect and the mixing to other channels. 

We can justify this procedure by evaluating the sus- 
ceptibility in an independent (and equivalent) method. 
By definition, in general, the susceptibility is to be iden- 
tified as the inverse of the potential curvature. For the 
purpose to compute the curvature inverse, we should con- 
sider the curvature matrix C whose dimensionless com- 
ponents are given by C uu = T 2 d 2 £lpN.ih/d(uu) 2 , C us = 
T 2 d 2 n PN j L /d(uu)d(ss),_C ue = T- l d 2 ^ mh /d(uu)d£, 
Cfi — T~ 4 d 2 ripNjh/d£d£, and so on. In the present case 
C is a 4 x 4 matrix. Then the diagonal components of 
C _1 give the susceptibility which is an involved expres- 
sion in terms of C uu , C us , C u g, etc. Roughly speaking, 
the diagonal part, C~^ 7 C" 1 , C^ 1 , corresponds to soft- 
mode propagators and the off-diagonal part, C us , C u i, 
Cu, Cn, and so on, corresponds to mixing vertices. It is 
immediate to make sure that C~ 1 certainly leads to ex- 
actly the same results as obtained from Eqs. (|20[) . (f2"Tj) . 
and (|22p . This matrix method has an advantage in giving 
us the mixing angle between each mode. 

As we can notice from Fig. [4] showing the susceptibility 
as a function of T, two crossovers associated with (uu) 
and I are located close to each other but do not coin- 
cide precisely. As long as we treat the chiral condensate 
and the Polyakov loop as independent variables as in the 
PNJL model, two crossovers attract each other to some 
extent but have a short "repulsion." Within this kind of 
model approach it is hence hard to explain the complete 
coincidence without fine tuning. 

One interesting strategy is not to explain the locking 
but to build a new model based on the complete locking 
of chiral restoration and deconfinement. As discussed 
in Ref . [l(J , most of lattice results support the idea that 
there is only one order parameter field cj> that is a mixture 
of the a meson and the electric glueball (Polyakov loop). 
Then, we could make a model with the chiral condensate 
given by (uu) oc 4>cos9 and the Polyakov loop by £ oc 
</>sin# with some potential energy for the mixing angle 
9 between them. The work along this direction is under 
progress [43j . 



D. Quark Number Susceptibility 

It is difficult to probe physical observables sensitive 
to the chiral and Polyakov loop susceptibility directly in 
experiments. In fact, it is impossible to count the num- 
ber fluctuation of the a meson and the glueball which 
eventually decay to the lightest n meson. From the ex- 
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T [MeV] 

FIG. 5: Quark number susceptibility calculated in the PNJL 
model for 2 + 1 flavor quark matter at zero density. 
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FIG. 6: Ratio of the fourth derivative to the quark number 
susceptibility calculated in the PNJL model for 2 + 1 flavor 
quark matter at zero density. 



perimental point of view the quark number susceptibility 
should be a better measure because the quark number is 
a conserved quantity. The fluctuation in the baryon mul- 
tiplicity would be directly related to the quark number 
susceptibility, Xg [3 Hi, EH- Also in Refs - EI SI Xg 
has been evaluated and discussed in the two-flavor PNJL 
model. Actually the PNJL model can reproduce Xq mea- 
sured on the lattice in the two-flavor case as beautifully 
illustrated in Ref. [ilj . 

We plot our results in the 2 + 1 flavor case in Fig. [5] 
We can see, as expected, that the 2+1 flavor quark 
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matter yields \ q greater than the two-flavor ease shown 
in Ref. [4l[. In the chiral limit Xq would scale as Nj, 
and thus the three-flavor value should be 3 2 /2 2 = 2.25 
times greater than the two-flavor value. Because s-quarks 
are massive in reality, this scale factor should become 
smaller. Let us choose one temperature to take an ex- 
ample for comparison. At the temperature T = 1.5T C ~ 
300 MeV, Fig. [5] reads around 2.7, while the two-flavor 
value is around 1.5, which leads to the ratio 1.8. This 
seems to be a reasonable number. 

It is interesting to define the following quantity (46[; 



Xi 4) = T 



^ 4 (-»PN.i L /r) 



(24) 



and take the ratio to the quark number susceptibil- 
ity. This ratio, /Xqi counts the number squared of 
quark content inside thermally excited particles carry- 
ing baryon number. Therefore, if quarks are liberated 
in the high temperature region, Xq /Xq — I 2 should fol- 
low, whereas the low temperature side should results in 
Xq^ IXq — 3 2 because of confinement. This is actually 
a general feature in the quasi-particle picture and at- 
tributed to the Boltzmann factor in the free fermionic 
partition function. 

Figure [5] shows this susceptibility ratio obtained in the 
2+1 flavor PNJL model. We see that the behavior per- 
fectly fits what is expected. A short conclusion that we 
should learn from this analysis is that X^q* IXq signifies 
the quark number but does not tell us whether the ther- 
mally excited particle is a confined nucleon or a set of 
three quarks. The latter is the case in the PNJL model. 



E. More on Thermodynamics 

Before proceeding into the finite density inquiry, we 
shall exemplify the success of the PNJL model by two 
more thermodynamic quantities. 

The trace of the energy momentum tensor is vanish- 
ing at the classical level when theory has no mass scale. 
We know that QCD in the chiral limit is scale invariant, 
which means that the trace of the energy momentum ten- 
sor in massless QCD is zero unless quantum corrections 
are taken into account. The QCD scale Aqcd arises from 
the dimensional transmutation due to the trace anomaly 
at the quantum level. 

In thermodynamics the traceless of the energy momen- 
tum tensor without mass gap means e — 3p — where p 
is the pressure given by — ^pnjl and e is the internal 
energy given by T 2 d{p/T) / dT . It is straightforward to 
evaluate e — Zp or the "interaction measure" from f2pN,iL 
in our model study. 

We show the model results in Fig. [71 The gross struc- 
ture with a peak around T c is in nice agreement with 
the recent lattice data (see Fig. 4 in Ref. [l!|). The 
peak height in the interaction measure is not as large 
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FIG. 7: Plot for the "interaction measure," i.e. (e — 3p)/T 4 , 
and the fermion contribution, Op M = [2m u d({uu) — {uu)o) + 
m a ((ss) - (ss)o)]. 
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FIG. 8: Plot for p/e as a function of e 1/4 . 



as that in Ref. [l9(, which is partly because of the fi- 
nite number of N T on the lattice and partly because of 
the smaller fermion contribution, 0p M , in our calculation. 
We present the results for 0p M = 2m u d((uu) — (uu)o) + 
m 3 ((ss) — (ss)o) also in Fig. [Jj We see that our results 
are significantly smaller than the results shown in Fig. 5 
in Ref. [19|. This is because the quark mass is different; 
the 7T mass in Ref. [l9[ is still around 220 MeV, while we 
choose m u d to yield the realistic it mass. 

We should be aware that the interaction measure, (e — 
3p)/T 4 , has only little to do with the trace anomaly in 
the PNJL model study. We have model inputs with mass 
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dimension, that is, the cutoff A. (There are four more 
dimensional parameters, <?s, gjy, °j & n d b but they can be 
all dimensionless in unit of A.) 

As a matter of fact, the peak structure is rather generic 
regardless of any specific model. One can understand this 
from the thermodynamic relation, 

£~3P _ T d (V J_\ , 90 

T 4 -^tKt't^J- [ ' 

The right-hand side is the temperature derivative of p/T 4 
where p/T 4 naively counts the effective degrees of free- 
dom v as plotted in Figs.[T]and[3] Therefore, so-called the 
trace anomaly, e — 3p, signifies how quickly v grows up as 
T increases. The pseudo-critical temperature is, by defi- 
nition, where v starts getting larger, and eventually v is 
saturated to the total number of particle species at high 
temperature. As a result the peak shape is inevitable 
associated with crossover behavior. It is not quite sur- 
prising in this sense that the PNJL model can mimic the 
trace anomaly in hot QCD around T c . 

In other words, it is the relation between e — 3p and the 
gluon condensate that is a non-trivial consequence from 
the trace anomaly. The interaction measure, (e — 3p)/T 4 , 
is governed not by the anomaly but by the thermodynam- 
ics which determines the gluon condensate in turn. 

Now that we have come by the pressure and the in- 
ternal energy, we can infer the sound velocity. Although 
the velocity of sound is given by c 2 — dp/de, the ratio 
p/e can approximate it in the high temperature limit. To 
compare our results to the available lattice data, we plot 
p/e as a function of e 1 / 4 in Fig. [5J which agrees quite 
well with Fig. 9 in Ref. [19(. We remark that the sound 
velocity has been investigated by means of the two-flavor 
PNJL model also in Ref. [13] where both of c 2 and p/e 
are presented. 



V. FINITE DENSITY RESULTS 

By adding one more axis in the direction of quark 
chemical potential we can investigate the order parame- 
ter behavior and the phase structure in wider perspective. 
In this work we limit ourselves to the chiral and decon- 
finement phase transitions and do not take account of the 
diquark condensation that plays an essential role in the 
color superconducting phase [4(| [H, IH, H3| ■ 




T [MeV] 



FIG. 9: Normalized light-quark chiral condensate {uu)/{uu)o 
in the fi-T plane, where (uu)o is the chiral condensate at 
T = n = Q. 

We remark that in the present parameter set the con- 
stituent quark masses turn out to be 

M ud = 336 MeV , M s = 528 MeV , (26) 

and the first-order phase transition is located at \i = 
345 MeV when T = Q, which is slightly above the light 
quark mass. This is a general feature to be explained in- 
tuitively. First, let us focus on the region at [i < 345 MeV 
where the system does not have any discontinuous tran- 
sition along the T — density axis. We can still locate 
the point where a non-vanishing baryon density appears 
at fi = M uc i, that is a sort of continuous phase transi- 
tion from the empty vacuum to degenerate quark matter. 
Next, once quark matter is concerned at [i ~ 345 MeV, 
the pressure of cold quark matter at a fixed value of 
fi becomes smaller for larger quark mass; for instance 
p oc /i 4 for massless quarks and p oc (/i 2 — M 2 ) 2 for mas- 
sive quarks. Thus, the kinetic energy favors lighter quark 
matter, that is, the chiral symmetric phase. The conden- 
sation energy gives a negative contribution to the pres- 
sure, that means that the chiral symmetric phase where 
the condensation energy is smaller is energetically favor- 
able again. In this way, one can expect that, as soon 
as the quark number density becomes substantial with fi 
going above M u d, the system tends to undergo a phase 
transition to the chiral symmetric phase. 



A. Chiral Phase Transition 

Figure[H]is a 3D plot for (uu) as a function of T and /i. 
We see that there is a discontinuity in the low tempera- 
ture and high density region, while the high temperature 
and low density region has continuous crossover. There- 
fore the phase diagram has an end-point of the first-order 
phase boundary, that is called the QCD critical point. 



B. Polyakov Loop 

It is interesting to see what happens in the Polyakov 
loop behavior on the fx-T plane. One may well anticipate 
that the coincidence of chiral restoration and deconfine- 
ment should persist in the finite density region. This 
expectation is partially true but too naive. We shall 
discuss the appropriate physical interpretation in what 
follows below. 
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FIG. 10: Polyakov loop I in the fi-T plane. 
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FIG. 11: Quark number density normalized by the free mass- 
less quark density, N c Nf(fi 3 /Sty 2 +T 2 /x/3), in the fi-T plane. 



We plot the Polyakov loop £ in the [i-T plane in Fig.fPJl 
It should be mentioned that we do not make another plot 
for £, for £ has a qualitatively same functional shape as I 
with small quantitative difference. 

From the comparison between the chiral condensate 
displayed in Fig. [5] and the Polyakov loop in Fig. [TDl we 
can readily perceive that two crossovers are linked in the 
entire region on the fi-T plane. For instance, we have 
already confirmed that two crossovers are simultaneous 
indeed at zero density in Fig. ® and we can find a first- 
order phase transition along the density axis at low tem- 
perature whose location is exactly the same in Figs.[5]and 
1101 The locking of chiral restoration and deconfinement 
remains at finite density in this sense. 

It would be misleading, however, to dive into a con- 
clusion that two phenomena of chiral restoration and de- 
confinement simultaneously take place in the high den- 
sity region. In view of the Polyakov loop behavior at 
low temperatures, in fact, the discontinuous jump is tiny 
and the expectation value of the Polyakov loop stays van- 
ishingly small even at fi > 345 MeV where chiral sym- 
metry is restored. Therefore, the discontinuous jump in 
the Polyakov loop signifies a first-order phase transition 
from nearly confined matter (£ ~ 0) with chiral symme- 
try breaking ((uu) ^ 0) to nearly confined matter {£ ~ 0) 
with chiral symmetry restoration ((iiii) ~ 0). 

It is an interesting question how the Polyakov loop be- 
haves such differently from the chiral condensates in the 
region of low temperature and high density. This is be- 
cause center symmetry is not broken at zero temperature 
even in the presence of dynamical quarks, and therefore, 
the expectation value of the Polyakov loop must stay van- 
ishing. The reason for preserved center symmetry is to be 
understood intuitively as follows; when the quark density 
is specified by a certain chemical potential, each energy 
level is occupied by a quark up to the Fermi surface. Be- 
cause quarks have color degeneracy, red, green, and blue 
quarks always sit on the same energy level, which makes a 
color singlet. One can easily see this really happening in 



the model from the Dirac determinant given in Eqs. (JTTJ) 
and (fT2| . That is, when /i > e we only have the second 
term out of the whole particle contribution, 



1 + e 3\s-»\/T + 3ee \e-„\/T + 3 J e a| 6 -„|/T > 



(27) 



that is exponentially dominant for large |e — [i\jT . This 
second term, e 3 ' e_AI '/ T , actually represents the three- 
quark occupation which does not couple £ and thus not 
break center symmetry. The third and fourth terms 
are one-quark and two-quark (which is equivalent to one 
anti-quark in color) contributions with non-trivial (non- 
singlet) color. Consequently, the Polyakov loop is insensi- 
tive to whether the quark degrees of freedom are present 
in the system or not. Strictly speaking, the PNJL model 
cannot deal with confinement, namely, nucleon wavefunc- 
tions as a bound state out of three quarks, but still, the 
low temperature region always has a signature of confine- 
ment (£ ~ 0) with respect to quarks. We would stress 
that this quarky confined phase is not a model artifact 
but physical one. We propose that this phase should be 
identified as the quarkyonic phase discussed in Ref. |5l|. 

One important suggestion emphasized in Ref. [5l| is 
that the baryon or quark number density serves as an 
order parameter to tell the quarkyonic phase from the 
hadronic phase. We have then calculated the quark num- 
ber density in the /.i-T plane to make a plot of Fig. QT] 
The quark number density is readily available from 



PNJL 



(28) 



In order to visualize in a sensible manner on the 3D plot, 
we normalize n q by the free massless quark density given 
by N c N f (fi 3 /37r 2 +T 2 ju/3), that is, in Fig. CCD we plot, 



AT = -2 

9 3/z 3 /^ 2 +3TV 

which should take a value from zero to unity. 



(29) 
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We can clearly confirm that the quark degrees of free- 
dom are relevant (i.e. N q ~ 0(1)) even in the region at 
r ~ and /i > M u d where £ ~ 0. This gives another 
evidence for our identification to the quarkyonic phase. 

It is a non-trivial finding from the present PNJL model 
study that N q surely behaves as an order parameter and 
locates the crossover point that coincides the chiral phase 
transition point. This coincidence is apparent at a glance 
of Figs. M and Q31 



C. Phase Diagram 

We now explore the phase diagram in the ji-T plane 
by the cross-section of Figs. [5] [TU1 and [11] at a certain 
height in the vertical axis. As we discussed, in the high 
density region in particular, the susceptibility peak po- 
sition does not make much sense, but the magnitude of 
the order parameter is a more suitable quantity to probe 
the physical state of matter. For example, the Polyakov 
loop susceptibility diverges at the critical end-point as 
well as the chiral susceptibility, but it does not result 
from the deconfinement transition but from the mixing 
to the chiral dynamics. Therefore we define the pseudo- 
critical temperature for w-quark chiral restoration by the 
condition: 



(uu) 



Wo 



T=T„Qx) 



(30) 



In the same way we can define the pseudo-critical tem- 
perature for s-quark chiral restoration by 



(ss) 



(ss) 



T—T s (fi) 



(31) 



Also, we can define the pseudo-critical temperature for 
deconfinement by 



T=T t (n) 



(32) 



Then, we can draw three distinct curves by T = T u (fi), 
T a (fi), Tt{n) on the fi-T phase diagram. The PNJL model 
prediction is shown in Fig. 1121 The solid curve represents 
T = T u (n) which is crossover in the low density region 
and turns a first-order phase transition in the high den- 
sity region accompanied by a critical end-point. We note 
that N q is nearly zero inside this solid curve. Because of 
explicit symmetry breaking by m s ^ 0, the T = T s (/j,) 
boundary is located at higher T and /i shown by the dot- 
ted curve in Fig. [T3J Of course, strictly speaking, chiral 
symmetry or even SUv(3) symmetry is not restored at 
any temperature or density, but (ss) can decrease up to a 
half of (ss)o smoothly. Actually this boundary hits T = 
and \x = 512 MeV which is not far from the constituent 
s-quark mass. In any case, the boundary by the dotted 
curve does not have a strong meaning because the change 
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FIG. 12: Phase diagram characterized by three quantities, 
namely, the u-quark chiral condensate, the s-quark chiral con- 
densate, and the Polyakov loop. 

in (ss) as a function of T is only gradual. Nevertheless, 
the region bounded by T u (fi) < T < T s (/i) is interesting. 
This is because the SUv(3) symmetry breaking becomes 
enhanced further in this region by chiral restoration for 
it-quarks and d-quarks but not for s-quarks (52J. 

It is surprising that the deconfinement crossover de- 
fined by the condition (|3"2"]) goes away from the chiral 
phase transition at finite density as indicated by the 
dashed curve which represents the T = Tg(fi) curve. As 
we have discussed, the Polyakov loop expectation value 
is always vanishing at zero temperature, and thus, this 
T = Ti(p) curve never hits the horizontal axis at T = 0. 
The region surrounded by T u (fi) < T < T^(/x) is what 
should be called the quarkyonic phase embodied in the 
PNJL model. 

As a final remark in this section we refer to the similar 
results presented in Figs. 16 and 17 in Ref. (4f| and the 
similar physical picture to the quarkyonic phase discussed 
in a different context in Ref. 15311. 



VI. QCD CRITICAL END-POINT 

We have already mentioned on the QCD critical end- 
point in the explanation of Fig. [12] The rest of this paper 
will be devoted to physics related to the QCD critical 
point. First of all, it is instructive to elucidate how the 
location of the critical point moves by the effect of the 
Polyakov loop. In the three-flavor NJL model with the 
Hatsuda-Kunihiro parameters, the location of the critical 
point is found to be 



(T E , ii e ) = (48 MeV, 324 MeV) , 



(33) 



in the three-flavor NJL model. The location is almost 
the same as in the two-flavor case. [See Ref. [55J for a 
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Chiral Suscept. 




FIG. 13: 3D plot for the light-quark chiral susceptibility Xud 
multiplied by (T/A) 2 in the jit-T plane. 

summary table and also Ref. The model depen- 

dence is nicely compiled also on Fig. 4 in Ref. [5q . In 
the three-flavor PNJL model the location goes up along 
the temperature to 

(T E , he) = (102 MeV, 313 MeV) , (34) 

in the present parameter set. This value is close to the 
two- flavor PNJL location first reported in my paper [l3| . 
The reason why the critical point moves toward higher 
temperature is that the artificial quark excitation at finite 
temperature and density is suppressed by the Polyakov 
loop average as exhibited in Fig. [3] and also discussed 
around Eq. (127]) . 

The question is to what extent we can trust the model 
prediction for the location of the QCD critical point or 
even its existence. In what follows we will discuss the 
dependence on ambiguous model parameters. So far, it 
is quite difficult to make any robust statement about the 
QCD critical point from model studies, that is our short 
conclusion. 



A. Divergent Susceptibility 

Before addressing the theoretical uncertainty on the 
QCD critical end-point, we will begin with standard ar- 
guments, that is, physical implication from the assump- 
tion that the QCD phase diagram holds a critical end- 
point. 

The importance of the QCD critical point lies in 
the fact that it is an exact second-order phase transi- 
tion point. Therefore the susceptibility diverges right 
at the end-point. Originally divergent growth in the 
chiral susceptibility Xu has been paid attention [571 ] 
which might lead to furious fluctuations in the a chan- 
nel and thus 7r fluctuations through the a «-» 2tt cou- 
pling. We have made a 3D plot in Fig. [13] to show the 
it-quark chiral susceptibility multiplied by (T/A) 2 , i.e. 
— jA~ 2 d 2 QpNjh/dm 2 id in the \i-T plane. We notice that 
the susceptibility has a singularity at the critical point. 




FIG. 14: 3D plot for the quark number susceptibility Xq mul- 
tiplied by (T/A) 2 in the \i-T plane. 



A physical quantity of more experimental interest is 
the quark number susceptibility. We plot Xq multiplied 
by (T/A) 2 in Fig. QH that is, -A~ 2 9 2 f7 PNJL /9/i 2 in the 
/i-T plane. Figure [14] shows a singularity at the QCD 
critical point which should translate into event-by-event 
fluctuations of baryon multiplicity. The global shape is 
just similar to that of the chiral susceptibility. It is, how- 
ever, different that the quark number susceptibility gets 
non-vanishing in the high temperature or density region 
whose behavior is closely linked to the quark number 
density in Fig. ITT] 



B. Columbia Diagram 

What we will reveal particularly in this work is the ro- 
bustness of the existence of the critical end-point, which 
is in part motivated by the lattice suggestion [58(. We 
can disclose another aspect of the phase diagram in the 
plane of the light and heavy quark masses [59( . Such a 
phase diagram is sometimes referred to as the "Columbia 
Diagram." 

The PNJL model results are summarized in Fig. [T5] 
Each curve represents the boundary between the first- 
order phase transition to the crossover. For m u d and m s 
below the curve, the chiral phase transition at finite T 
is of first-order, and otherwise, it is crossover. We draw 
a line m s /m u d = 24.67 which crosses the physical point, 
and add two lines at m u d = and m s = 0, respectively, 
for the eye guide. 

This plot poses us a problem in the model study based 
on the NJL-type description. It is that the first-order 
phase transition region at fi = is substantially smaller 
than what has been observed in the lattice QCD simu- 
lations. In the m s = case, for instance, the critical 
value of the light-quark mass is m u — 2.1 MeV in this 
work, and on the m u d = axis, the critical strange quark 
mass is m s = 8.8 MeV, which are smaller by one or- 
der of magnitude at least as compared to the lattice em- 
pirical value (58|. This fact implies that the first-order 
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phase transition with massless three flavors is presumably 
weaker in the NJL model than realistic. Then, the criti- 
cal end-point at physical quark mass cannot avoid being 
far away from zero density. That is, the density has to 
increase significantly until the boundary eventually hits 
the physical mass point. This is why the NJL- type model 
has common tendency to lead to the critical end-point at 
relatively high density above fi ~ 300 MeV. 

Because the PNJL model predicts the QCD critical 
point, the first-order transition region expands with in- 
creasing /i as shown in Fig. 1151 The boundary surface 
is thus standard but not quite consistent with the recent 
lattice observation [58| . This is problematic to the model 
study if the lattice results are correct. The model study 
has, however, unknown factors which could make a dras- 
tic change in the order. Here, we will point out two major 
effects; one is the Ua(1) anomaly reduction in a medium 
and the other is the induced vector-channel interaction. 



1. Anomaly strength 



FIG. 16: First-order transition boundary depending on the 
strength of the 't Hooft coupling go, where go is a value 
fixed in the vacuum in the Hatsuda-Kunihiro parameter set. 
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We have already noted that the first-order transition 
region on the Columbia diagram obtained in the PNJL 
model is significantly smaller than the lattice results. 
One likely explanation for this is that the 't Hooft (six- 
quark) coupling constant, <?d, is weaker in the NJL model 
estimate than realistic because of cutoff artifact. It 
should be mentioned that the value of gu is fixed to re- 
produce the r]' mass, which is as large as 957.5 MeV and 
is greater than the cutoff A = 631.4 MeV. It is not un- 
likely that the strength of c/d has been underestimated 
to reproduce such a large mass in this cutoff model. 

Figure [TBI shows the first-order transition boundary on 
the m u d-m s plane. Here ^do denotes the standard value 
in the Hatsuda-Kunihiro parameter set. Because <?do has 
been fixed in the vacuum, g^ in a hot and dense medium 
might take a different (presumably smaller) value. As 



FIG. 17: Dependence of the critical point location on the 
strength of the 't Hooft coupling constant. 



we anticipated, the first-order region becomes wider with 
larger g^ and narrower with smaller <?r> It should be 
noted that the plot is made in the linear scale in Fig. [TBI 
while the scale is logarithmic in Fig. [151 

One can then expect that the QCD critical point 
should move accordingly as gn changes. We show the 
location of the QCD critical point for various values of 
<7d in Fig- E3 We can learn two lessons from this fig- 
ure: One is that the QCD critical point can be located 
at higher temperature and lower density if g^> is under- 
estimated in the NJL model study due to too heavy rf 
out of model reliability. The other is in a sense oppo- 
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FIG. 18: Dependence of the critical point location on the 
strength of the vector-channel interaction. 



site to the first one. The QCD critical point might be 
absent from the QCD phase diagram if <7d is suppressed 
by the effective Ua(1) restoration at high density. Ac- 
tually, only 35% reduction is enough to make the QCD 
critical point disappear from the phase diagram. If the 
suppression is exponential like gn(fj) = e~^ 2 / fJ, og DQ [l^ |, 
35% reduction is within a reasonable reach. 

Then, one could change the scenario in Fig. [IS] by in- 
troducing a ^-dependent value for c/d- For instance, if 
one assumes an exponential ansatz, <7d(m) = °9t>qi 
one could find some fiQ that produces a boundary surface 
with bending behavior that the first-order transition re- 
gion shrinks with increasing /i. 



2. Vector- channel interaction 

It is not only the Ua(1) anomaly term but also the 
vector-channel interaction term in Eq. (Q that can affect 
the location of the QCD critical point. We remark that 
Cy does not break chiral symmetry at all, and besides, 
the zeroth component corresponds to the density oper- 
ator (ip^ip) 2 . Therefore, it is conceivable to expect that 
the finite density environment brings about non-zero gy 
even though we choose gy to be zero in the vacuum. 

There is no constraint at all for the choice of induced 
gy at finite density. We have no knowledge on even its 
sign. Since we regard gy in the present study as induced 
in dense quark matter, the choice of gy has nothing to 
do with the vector meson property in the vacuum. [It 
might be related to an in-medium modification.] It is 
presumably appropriate to measure the strength of gy in 
unit of gs, and we just try various values of gy to grasp 
a feeling of its effect. 



There are two modifications necessary to accommodate 
the vector-channel interaction. The condensation energy 
should be O con d — > f2 C ond— gyn q where we already defined 



in Eq. 



At the same time, the quark chemical 



potential should be replace by the renormalized one, 

H r = /i - 2gyn q , (35) 

like the quark mass replaced by the constituent one. 
Then, we have to solve the number constraint equation, 
n q = — <9f2pNjL/<9/i r , together with the four gap equa- 
tions self-consistcntly. In view of the condensation energy 
part, positive gy seems to decrease the free energy for 
non-zero n q , but the chemical potential renormalization 
overcomes it and the free energy becomes greater. Be- 
cause chiral symmetric phase has smaller M u d and thus 
larger n q , the vector-channel interaction with gy > de- 
lays chiral restoration. 

The results are summarized in Fig. llSl in the same way 
as in Fig. [T7] It is remarkable that the qualitative feature 
is quite similar to Fig. 1171 Thus, we can draw the same 
conclusions as in the case of the Ua(1) anomaly term. 

The QCD critical point could be absent again in the 
case when gy is greater than around 0.206gs- The critical 
value turns out to be consistent with the known value in 
the two-flavor case [HI, [33[ . It should be noted that the 
normalization of gs in Refs. [32l. |33| is different from the 
present convention by a factor 2. 

This value of critical gy is small as compared to the 
empirical value suggested by the Fierz transformation. If 
we take care of the effect of the effective Ua(1) restora- 
tion, as we illustrate in Fig. [T71 the critical gy could be 
even smaller. 



VII. CONCLUSIONS 

We have formulated the 2 + 1 flavor PNJL model with 
a simple ansatz for the Polyakov loop effective potential. 
We first confirmed that our model setup works pretty 
well to account for recent results in the zero-density lat- 
tice QCD simulation. We then explored our perspective 
toward the finite-density QCD phase transition. 

The phase diagram in our model study turns out to 
have three (crossover) boundaries corresponding to ud- 
quark chiral restoration, s-quark chiral restoration, and 
deconfincment characterized by the Polyakov loop expec- 
tation value. We have also computed the quark number 
density and found that its behavior is governed by the it- 
quark chiral condensate. Our phase diagram is consistent 
with the large N c argument and, in particular, we iden- 
tified the phase region with vanishing Polyakov loop and 
nonzero quark number density as the quarkyonic phase. 

It would be intriguing to include the diquark conden- 
sates to describe a family of the color superconducting 
phases. The large N c argument cannot access physics of 
color superconductivity, and thus, nothing so far could 
predict the fate of the quarkyonic phase region under 
the effect of color superconductivity. One possibility is 
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that the quark-hadron continuity realizes at low tem- 
perature and high density, and there appears crossover 
from the quarkyonic phase to the color superconducting 
phase, which is to be interpreted as crossover from con- 
fined to deconfincd quark matter. We remark that this 
continuity scenario requires three flavors and there have 
already been other scenarios (i.e. phase transitions with 
the coexisting regions) within the PNJL model frame- 
work [ii Hi|. 

Also, we have closely investigated parameter depen- 
dence of the location of the QCD critical point. We 
demonstrated that the QCD critical point moves quite 
easily in accord to the choice of the Ua(1) anomaly 
strength and the vector-channel interaction. Both are 
not under theoretical control at finite density. In fact, 
we have found that the critical values of these parame- 
ters are within a conceivable range in dense quark matter. 
That means, not only the location but also the existence 
of the QCD critical point is not robust at all in the model 
study. 

Although we limited our discussions only to numeri- 
cal results in this paper, it could be viable to examine 
the density dependence of the Columbia diagram in an 
analytical way [60]. Analytical understanding should be 
useful for the lattice QCD study from the zero density 
approaching toward the critical point. 

To establish the existence or non-existence of the QCD 



critical point, anyway, we must wait for future develop- 
ment of the finite-density lattice simulation, or experi- 
mental confirmation. [See also Refs. [fill, [62j for proposed 
experimental signatures.] 

Finally, let us refer to some of recent attempts in the 
PNJL model. The neutrality condition has been con- 
sidered in Ref. [63j |. It is known that the neutrality 
plays an important role especially in bulk superconduct- 
ing quark matter. Because the (untraced) Polyakov loop 
is a color matrix, there arise non-trivial coupling between 
the Polyakov loop and the color chemical potential, which 
may bring a subtle difficulty involving the gauge choice. 
This is a future problem. In Ref. [HJl^] the imaginary 
chemical potential has been considered. This is also an 
interesting future direction toward the nature of high- 
density quark matter. 
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